source('00_util_scripts/mod_seurat.R')
source('00_util_scripts/mod_bplot.R')

sobj <- read_rds('CRC-I/data/seekgene/crc-starsolo-annotated.rds')

sg_meta <- sobj@meta.data |>
  as_tibble(rownames = '.cell')

# find most frac changed cell type between II+IT vs TT --------
discov_frac_change2 <- function(x, y){
  x |>
    filter(zhang2020 == 'y') |>
    discov_frac_change(varoi = genotype, typeoi = latent_cluster,
                       var.1 = TT, var.2 = IT)
}

## myeloid ------------
frac_lfc_myl <- sg_meta |>
  discov_frac_change2('Myeloid cell')

frac_lfc_myl |>
  ggplot(aes(subtype, log2fc_frac, fill = type)) +
  geom_col() +
  coord_flip() +
  theme_pubr(legend = 'right') +
  labs_pubr() +
  scale_fill_manual(values = c('blue','red','grey')) +
  labs(title = 'Myeloid cell type fraction changes TT vs II+IT',
       y = 'log2fc of fraction change TT vs II+IT')

sobj_myl <- sobj |>
  filter(zhang2020_main == 'Myeloid cell')

Idents(sobj_myl) <- sobj_myl$latent_cluster

all_marker_myl <- sobj_myl |>
  FindAllMarkers()

myl_s100b_deg <- all_marker_myl |>
  as_tibble() |>
  filter(str_detect(cluster, 'S100B') & p_val_adj < .05)

myl_s100b_orago <- myl_s100b_deg |>
  filter(avg_log2FC > 0) |>
  pull(gene) |>
  enrichGO(OrgDb = 'org.Hs.eg.db',
           keyType = 'SYMBOL',
           ont = 'BP',
           readable = TRUE)

myl_s100b_orago@result |>
  head(10) |>
  plot_enrichment()

myl_s100b_orakegg <- myl_s100b_deg |>
  filter(avg_log2FC > 0) |>
  pull(gene) |>
  bitr(fromType = 'SYMBOL', toType = 'ENTREZID',
       OrgDb = 'org.Hs.eg.db') |>
  pull(ENTREZID) |>
  enrichKEGG()

myl_s100b_orakegg@result |>
  head(10) |>
  plot_enrichment()

## CD4 T -------------
sg_meta |> count(zhang2020_main)

frac_lfc_cd4 <- sg_meta |>
  discov_frac_change2('CD4 T cell')

frac_lfc_cd4 |>
  ggplot(aes(subtype, log2fc_frac, fill = type)) +
  geom_col() +
  coord_flip() +
  theme_pubr(legend = 'right') +
  labs_pubr() +
  scale_fill_manual(values = c('blue','red','grey')) +
  labs(title = 'CD4+ T cell type fraction changes TT vs II+IT',
       y = 'log2fc of fraction change TT vs II+IT')

## CD8 T -------------
frac_lfc_cd8 <- sg_meta |>
  discov_frac_change2('CD8 T cell')

frac_lfc_cd8 |>
  ggplot(aes(subtype, log2fc_frac, fill = type)) +
  geom_col() +
  coord_flip() +
  theme_pubr(legend = 'right') +
  labs_pubr() +
  scale_fill_manual(values = c('blue','red','grey')) +
  labs(title = 'CD8+ T cell type fraction changes TT vs II+IT',
       y = 'log2fc of fraction change TT vs II+IT')

## B cell -------------
frac_lfc_bc <- sg_meta |>
  discov_frac_change2('B cell')

frac_lfc_bc |>
  ggplot(aes(subtype, log2fc_frac, fill = type)) +
  geom_col() +
  coord_flip() +
  theme_pubr(legend = 'right') +
  labs_pubr() +
  scale_fill_manual(values = c('blue','red','grey')) +
  labs(title = 'B cell type fraction changes TT vs II+IT',
       y = 'log2fc of fraction change TT vs II+IT')

## crc4 --------
frac_lfc_crc4 <- crc4@meta.data |>
  mutate(istt = ifelse(genotype == 'TT', 'TT', 'other')) |>
  discov_frac_change(istt, manual_main, TT, other)
  
frac_lfc_crc4 |>
  ggplot(aes(subtype, log2fc_frac, fill = type)) +
  geom_col() +
  coord_flip() +
  theme_pubr(legend = 'right') +
  labs_pubr() +
  scale_fill_manual(values = c('blue','red','grey')) +
  labs(title = 'Immune cell type fraction changes TT vs II+IT',
       y = 'log2fc of fraction change TT vs II+IT')


# Augur to assess TT disturb what cell type most --------
skg_auc <- sobj |>
  filter(!is.na(latent_cluster)) |>
  calculate_auc(label_col = 'genotype',
                cell_type_col = 'zhang2020_main')

skg_auc |>
  plot_lollipop()

sj.main <- sobj |>
  filter(!is.na(latent_cluster)) |>
  SplitObject(split.by = 'zhang2020_main')

sub_auc <- sj.main |>
  map(\(x)calculate_auc(x, label_col = 'genotype',
                        cell_type_col = 'latent_cluster'),
      .progress = T)

sub_auc$`B cell` |>
  plot_lollipop() +
  labs(x = 'subtype', title = 'B cell') +
  theme_pubr()

sub_auc |>
  map('AUC') |>
  list_rbind(names_to = 'main') |>
  write_csv('CRC-I/results/seekgene.augur.csv')

skg_auc$AUC |>
  write_csv('CRC-I/results/seekgene.main.augur.csv')

## crc4 -------
crc4 <- read_rds('CRC-I/data/crc_merge4_immune.rds')

augur_crc4.tt <- crc4 |>
  mutate(tt = genotype == 'TT') |>
  calculate_auc(label_col = 'tt',
                cell_type_col = 'manual_main', n_threads = 7)

augur_crc4.tt |>
  plot_lollipop() +
  theme_pubr() +
  labs(x = 'Cell type')

## crc3 ---------
crc3 <- read_rds('CRC-I/data/crc_merge3.rds')

crc3.meta <- read_csv('CRC-I/results/crc_merge3_myl.namely.csv')

crc3.meta

crc3 |>
  filter()

# MIF expressing myeloid cell? ----
sobj_myld <- sobj |>
  filter(zhang2020_main == "Myeloid cell") |>
  GetAssay() |>
  CreateSeuratObject()

sobj_myld@meta.data <- sobj@meta.data |>
  filter(zhang2020_main == "Myeloid cell")

sobj_myld <- sobj_myld |>
  quick_process_seurat(pcs = 25, res = 2) |>
  FindClusters(algorithm = 4)

Idents(sobj_myld) <- 'latent_cluster'

mif_clus <- sobj_myld |>
  FindAllMarkers(features = 'LGALS9', only.pos = T) |>
  pull(cluster)

g1 <- sobj_myld |>
  DimPlot(group.by = 'latent_cluster', cols = DiscretePalette(25))

g2 <- sobj_myld |>
  mutate(galectin9_positive = latent_cluster %in% mif_clus) |>
  DimPlot(group.by = 'galectin9_positive',
          cols = c('grey','purple'))

g1 + g2

g3 <- frac_lfc_myl |>
  mutate(latent_cluster = as.character(subtype)) |>
  select(latent_cluster, type) |>
  left_join(sobj_myld, y = _) |>
  DimPlot(group.by = 'type',
          cols = c('blue','red','grey')) +
  ggtitle('Proportional changes')

g1 + g3

gal9hi_deg <- frac_lfc_myl |>
  mutate(latent_cluster = as.character(subtype)) |>
  select(latent_cluster, type) |>
  left_join(sobj_myld, y = _) |>
  FindMarkers(ident.1 = 'Increased in TT', group.by = 'type') |>
  as_tibble(rownames = 'gene')

gal9hi_deg |>
  filter(avg_log2FC < .05) |>
  write_csv('CRC-I/results/galectin9_mf_deg.csv')

library(ggrepel)

gal9hi_deg |>
  plot_bill_volc('Galectin-9+ macrophage')

gal9_upgo <- gal9hi_deg |>
  filter(p_val_adj < .05 & avg_log2FC > 0) |>
  pull(gene) |>
  enrichGO(OrgDb = 'org.Hs.eg.db', keyType = 'SYMBOL', ont = 'BP')

gal9_upgo <- gal9_upgo |>
  clusterProfiler::simplify()

gal9_upgo@result |>
  plot_enrichment(n = 9)

gal9_upkegg <- gal9hi_deg |>
  filter(p_val_adj < .05 & avg_log2FC > 1) |>
  pull(gene) |>
  bitr(fromType = 'SYMBOL', toType = 'ENTREZID', OrgDb = 'org.Hs.eg.db') |>
  pull(ENTREZID) |>
  enrichKEGG()

gal9_upkegg@result |>
  plot_enrichment()

# TRPM2-hi macrophage -----------
sobj <- read_rds('CRC-I/data/crc_merge4_immune.rds')

sobj$hpca_main_c4 |> unique()

sobj |> DotPlot('TRPM2', group.by = 'zhang.main', cols = 'RdYlBu',
                cluster.idents = T)

m2.fine.expr <- last_plot() |>
  pluck('data')

m2.fine.expr |>
  mutate(id = fct_reorder(id, avg.exp)) |>
  ggplot(aes(features.plot, id, color = avg.exp.scaled, size = pct.exp)) +
  geom_point() +
  theme_bw() +
  labs(x = 'Gene', y = 'Cell type', title = 'TRPM2 in CRC') +
  scale_color_gradient2(low = 'blue', high = 'red') +
  labs_pubr()

sobj.myl <- sobj |> filter(manual.main == 'Myeloid_cells')

sobj.myl$seurat_clusters |> unique()

sobj.myl |> DotPlot('TRPM2', cols = c('RdYlBu'))

sobj.myl |> ggplot(aes(seurat_clusters, fill = zhang.fine)) +
  geom_bar()

sobj.myl <- sobj.myl |>
  quick_process_seurat(c('orig.ident', 'dataset'), skip_norm = T)

zh20.ref <- read_rds('CRC-I/data/crc10x_singler_ref.rds')

sobj.myl <- sobj.myl |>
  mark_cell_type_singler(zh20.ref, fine_label = T,
                         sc_ref = T, new_label = 'zhang.fine2')

sobj.myl |>
  mark_cell_type_singler(zh20.ref, fine_label = T,
                         sc_ref = F, new_label = 'zhang.fine2')

sobj.myl |> DimPlot(group.by = 'zhang.fine2', cols = 'Paired')

sobj.myl |> DotPlot('TRPM2', cols = c('RdYlBu'), group.by = 'zhang.fine2')

sobj.myl <- sobj.myl |>
  mutate(manual.fine = case_when(str_detect(zhang.fine, 'pDC') ~ 'pDC',
                                 str_detect(zhang.fine, 'cDC') ~ 'cDC',
                                 str_detect(zhang.fine, 'Mono-') ~ 'Mono',
                                 .default = 'Macrophage'))
